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The large systems of complex linear equations that are generated in QCD problems often have 
multiple right-hand sides (for multiple sources) and multiple shifts (for multiple masses). Deflated 
GMRES methods have previously been developed for solving multiple right-hand sides. Eigen- 
vectors are generated during solution of the first right-hand side and used to speed up convergence 
for the other right-hand sides. Here we discuss deflating non-restarted methods such as BiCGStab. 
For effective deflation, both left and right eigenvectors are needed. Fortunately, with the Wilson 
matrix, left eigenvectors can be derived from the right eigenvectors. We demonstrate for difficult 
problems with kappa near k c that deflating eigenvalues can significantly improve BiCGStab. We 
also will look at improving solution of twisted mass problems with multiple shifts. Projecting 
over previous solutions is an easy way to reduce the work needed. 
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1. Introduction 

Current lattice QCD simulations attempt to reach physical up and down quark masses. In 
this regime, standard linear solvers used in quark propagator calculations converge very slowly. 
Roughly speaking, the rate of convergence is proportional to the square root of the ratio of the 
smallest eigenvalue to the largest eigenvalue of the Dirac matrix. A remedy of this problem is to 
deflate some of the eigenvectors corresponding to the smallest eigenvalues |JlJ]. For restarted GM- 
RES, this was done by augmenting the Krylov subspace with approximate eigenvectors with small 
eigenvalues. The resulting algorithm is called GMRES with deflated restarting or GMRES-DR [0]. 
One advantage of GMRES-DR is that eigenvectors are calculated simultaneously while solving the 
linear system and no separate calculation is needed. The eigenvectors calculated are approximate 
and their accuracy increases with each restart. In addition, eigenvectors computed with GMRES- 
DR could be used to accelerate the convergence for subsequent right-hand sides. This is a common 
situation in lattice QCD calculations where one needs to find the quark propagator from all lat- 
tice sites using noise methods. The algorithm is called GMRES -Proj and is based on combining 
restarted GMRES with a projection over previously determined eigenvectors [^J. For multiple right 
hand sides the following two main steps are used: 

• Solve the first right hand side using GMRES-DR. 

• For subsequent right-hand sides, solve by alternating between a minimal residual projection 
step over the right eigenvectors with smallest eigenvalues obtained from GMRES-DR and 
one or more cycles of GMRES. 

Here, we extend this approach for multiple right-hand sides by replacing GMRES in the second 
step with BiCGStab. Since BiCGStab is a non-restarted method, the deflation step will be applied 
only once in the beginning. It will be used to obtain a better initial guess for the solution. In 
addition, with right eigenvectors the Min.Res. projection does not do a good enough job of reducing 
the crucial components in the direction of eigenvectors corresponding to small eigenvalues. So, 
instead, we will project using both right and left eigenvectors. In the case of Wilson fermions, left 
eigenvectors are obtained from right eigenvectors because the relation 75^^/5 = is satisfied by 
the Wilson Dirac operator Ayr- 

2. Algorithm for D-BiCGStab(k) 

Assuming that the first right-hand side was solved using GMRES-DR giving both the solution 
as well as approximate right eigenvectors. The deflated BiCGStab for subsequent right-hand sides 
given in [||] is as follows: 

• Consider the system Ax = b. Let xo be an initial guess and ro = b —Axo be the initial residual. 

• Let vi,V2,...,V£ an orthonormal basis for the set of k right eigenvectors and u\,U2,...,Uk an 
orfhonormal basis for the k left eigenvectors. Define U as n x k matrix whose columns are 
the left eigenvectors U = [ki,W2, .. Similarly, V is n x k matrix whose columns are the 
right eigenvectors V = [v\,Vz, ...,vjt]. 
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• Solve the k x k linear system U^AVy = Wr$ for y and construct an improved initial guess 
4 ew =x + Vy. 

• Apply BiCGStab to solve the system using xff w as initial guess. 
3. Results for D-BiCGStab(k) 

Deflated BiCGStab is first tested on quenched configurations generated using the Wilson plaquette 
action at /3 = 6.0 on 16 and 20 3 x 32 lattices (see §\§ for results on dynamical configurations). 
For Wilson fermions, we tune fc to be close to the critical value K c in order to make it a difficult 
but physical problem. The value of Kc was determined on each configuration from the condition 
that the real part of the smallest eigenvalue of A vanishes. In the following we solve the even-odd 
preconditioned system.The first right-hand side is solved using GMRES-DR(m,k) where k is the 
number of deflated eigenvectors and m is the maximal dimension of the subspace. Both m and k are 
varied but the difference m — k = 20 is kept fixed. The accuracy of the eigenvectors generated with 
GMRES-DR is increased by reducing the value of the residual norm at convergence when solving 
the first right-hand side. In the following, the first right hand side was solved with ratio of the norm 
of the residual at convergence to the norm of the initial residual 10~ 8 , 10~ 10 and 10~ 14 . This will 
be denoted by LI, L2 and L3 respectively with L3 corresponding to most accurate eigenvectors. 
For the second right-hand side, the relative residual norm at convergence is required to be 10~ 8 . 
For comparison, the second right-hand side is solved using standard BiCGStab and GMRES(20)- 
Proj(k). When convergence to the desired relative residual norm is not reached, this is indicated by 
the letter "F" in tables. Results are shown in Tables |l|,^|for a sample of three configurations. 



c# 




m,k 


GMRES -DR(m,k) 


D-BiCGStab(k) 


BiCGStab 








I st rhs 


LI 


L2 


L3 




1 


0.158383 


30,10 


650 


746 


708 


490 


720 






40,20 


560 


726 


550 


472 








50,30 


550 


764 


590 


402 








60,40 


540 


692 


504 


356 




2 


0.158399 


30,10 


710 


818 


718 


644 


806 






40,20 


620 


832 


584 


486 








50,30 


610 


796 


542 


370 








60,40 


600 


774 


458 


418 




3 


0.157924 


30,10 


710 


424 


420 


412 


776 






40,20 


620 


382 


388 


308 








50,30 


610 


500 


356 


306 








60,40 


600 


398 


374 


270 





Table 1: Results for the 16 4 lattice. Matrix-vector products listed correspond to relative residual norm 10 8 . 

Comparing the results for D-BiCGStab(k) to BiCGStab, we find that the deflation step leads 
to a considerable improvement and occasionally to a "breakthrough" as with the second configura- 
tion of the 20 3 x 32 lattice (Table ^| with (m,k) = (70,50) and accurate eigenvectors). For a fixed 
number of deflated eigenvectors, it is found that the improvement increases as the accuracy of the 
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c# 




m,k 


GMRES -DR(m,k) 


D-BiCGStab(k) 


BiCGStab 








V'rhs 


LI 


L2 


L3 




1 


0.157200 


30,10 


2310 


1974 


1738 


1244 


2332 






40,20 


1980 


2396 


1532 


1170 








50,30 


1590 


988 


788 


830 








60,40 


1540 


1414 


844 


646 








70,50 


1510 


978 


842 


676 




2 


0.157044 


30,10 


F 


1690 


1690 


1690 


1812 






40,20 


F 


1700 


1700 


1700 








50,30 


1410 


794 


812 


694 








60,40 


1360 


878 


720 


584 








70,50 


1330 


850 


748 


282 




3 


0.157095 


30,10 


1710 


1524 


1214 


1136 


1534 






40,20 


1220 


1188 


1134 


912 








50,30 


1170 


1172 


940 


918 








60,40 


1120 


1150 


912 


808 








70,50 


1110 


1214 


842 


656 





Table 2: Results for the 20 3 x 32 lattice. Matrix-vector products listed correspond to relative residual norm 

io- 8 . 



eigenvectors increases. So, the least number of matrix-vector products will correspond to the L3 
columns. As expected also, the more eigenvectors we deflate the smaller the number of matrix- 
vector products. However, after certain optimal number of eigenvectors, an increase in the number 
of the deflated eigenvectors does not lead to large improvement. This optimal number of eigenvec- 
tors was found to increase as the volume of the lattice increases but fortunately not linearly. Note 
also that when GMRES-DR does not converge, we don't get good eigenvectors and little improve- 
ment is obtained by deflation (see results for the second configuration in Table ||). For illustration, 
results for the first configuration of the 20 3 x 32 lattice are shown in Figs. [l|-||. In Fig. [I], a closer 
look at the eigenvalue spectrum near the origin shows a very small eigenvalue as well as other small 
eigenvalues. In Fig. ^, we show the effect of increasing the number of deflated eigenvectors as well 
as the effect of using more accurate eigenvectors. In Fig. we compare results for BiCGStab, 
GMRES-DR(50,30), GMRES(20)-Proj(30) and D-BiCGStab(30) using the most accurate vectors. 
Both deflated BiCGStab and GMRES -Proj improve substantially over BiCGStab. 

4. A note on multi-mass solvers for Twisted-Mass QCD 

In Twisted Mass QCD, quarks are introduced as pairs with a modified mass term [01]. The 
fermionic part of the action for a degenerate pair of quarks is given by (see [|J] for unexplained 
notations): 

5 / = LL^{i[rv( V v + V v)-V;V v ]+m + / J U7 5 T 3 }^, (4.1) 
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Figure 1: Small eigenvalues near the origin for the first configuration of the 20 x 32 lattice. 
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Figure 2: 

deflated ei 



Results for quenched Wilson configuration on 20 3 x 32. Left: effect of increasing the number of 
genvectors with L3 accuracy. Right: Effect of increasing the accuracy of the deflated eigenvectors. 



where *F = (u , d) is a doublet of the up and down quarks and m and pL are the standard and 
twisted mass terms respectively. Because of the twisted mass term, it is not possible to apply multi- 
mass solvers to twisted mass problems simultaneously with even-odd preconditioning. Multi-mass 
solvers as CGS can be used if even-odd preconditioning is not implemented [^]. In the following we 
describe how one can accelerate the convergence of twisted-mass problems with multiple masses 
and even-odd preconditioning. The method is based on solving the systems serially but using an 
improved initial guess by making a minimal residual projection over available solutions of the 
previous systems. It is described as follows: 

• Consider the systems Apt; = b where A, is the even-odd preconditioned Twisted-Mass Dirac 
operator that corresponds to K:f,/i,. We assume the system at maximal twist, so fcf is the 
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Figure 3: Comparing results for the first configuration of the 20 3 x 32 lattice 

critcal K value that corresponds to the twisted mass value /I,-. 

• Assume that we solved the systems for i = l,..k. In order to accelerate the solution of the 
system k + 1, we perform a minimal residual projection over the k previous solutions as 
follows: 

- let x Q k+l be the initial guess and r° +1 = b — Ak+ix k+l the initial resiual. Define Q as the 
n x k matrix whose columns are the previous k solutions, i.e. Q = [xi ,X2, ■■■,Xk]- 

- Solve the k x k system A^^A^iQd = G t ^J +1 ^ +1 for d and construct an improved 
initial guess = x® +l + Qd. 

The method was tested on 20 3 x 32 lattice with quenched configurations at /3 = 6.0 with the Wilson 
plaquette action. The twisted mass fermion action for a degenerate doublet of quarks at maximal 
twist for 1 1 values of /J. is used. The values of K c for each value of /I is determined using a linear 
fit of the four (K" f ,/x) pairs in Three parameters affect the performance of the method. First 
is the separation between successive masses, second the number of available solutions to project 
over and third whether we solve the heaviest or the lightest mass first. In Table |[ we compare the 
number of matrix-vector products when zero initial guess is used to the case where projection over 
previous solutions is done for a typical mass separation. Although the projection step is done only 
once per shift, we found a considerable reduction in the total number of matrix-vector products. 
The overall reduction seems to be similar whether we solve the heaviest or the lightest mass first 
for that mass separation. 

5. Conclusions 

For problems with multiple right-hand sides, a combination of GMRES-DR for the first right- 
hand side and a deflated BiCGStab for subsequent right-hand sides was tested on typical lattice 
volumes. It was found to give a considerable reduction of the matrix-vector products by a fac- 
tor of approximately 5 for the L3 case. The improvement level increases as the accuracy of the 
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Mass Number 






x° = 


With projection 
high — > low 


With projection 
low^high 


1 


0.157290 


0.005 


1270 


1000 


1270 


2 


0.157210 


0.009 


1150 


730 


1030 


3 


0.157130 


0.013 


1030 


550 


880 


4 


0.157050 


0.017 


880 


430 


700 


5 


0.156970 


0.021 


790 


400 


520 


6 


0.156890 


0.025 


730 


280 


400 


7 


0.156810 


0.029 


640 


280 


310 


8 


0.156730 


0.033 


580 


310 


220 


9 


0.156650 


0.037 


520 


340 


190 


10 


0.156570 


0.041 


490 


400 


160 


11 


0.156490 


0.045 


460 


460 


190 


Total MVP 






8,540 


5,180 


5,870 



Table 3: Effect of the projection step with GMRES-DR(40,10) on 20 3 x 32 lattice starting from highest to 
lowest mass or vice versa with Aji = 0.004. 



eigenvectors increase. Deflated BiCGStab was tested with left-right projection. The method does 
not add extra work for the case of Wilson fermions where the left eigenvectors are related to the 
right eigenvectors through 75 multiplication. For Twisted-Mass problems with multiple shifts and 
even-odd preconditioning, it was found that improving the initial guess using a minimal residual 
projection over previous solutions reduces the matrix- vector products by a factor of 20% — 50%. 
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